In [1]:
import numpy as np
import pandas as pd
import scanpy.api as sc
sc.settings.verbosity = 3 # verbosity = 3: errors, warnings, info, hints
sc.settings.set_figure_params(dpi=70) # dots (pixels) per inch determine size of inline figures
sc.logging.print_version_and_date()
from matplotlib import pyplot as plt
import matplotlib.colors as mcol
cm1 = mcol.LinearSegmentedColormap.from_list("BuRd",["b","r"])
Loading Data¶
In [2]:
data_folder = "data/regression/"
filename_data = data_folder+"matrix.mtx"
filename_genes = data_folder+"genes.tsv"
filename_barcodes = data_folder+"barcodes.tsv"
adata = sc.read(filename_data, cache=True).transpose()
adata.var_names = np.genfromtxt(filename_genes, dtype=str)[:, 1]
adata.smp_names = np.genfromtxt(filename_barcodes, dtype=str)
Preprocessing¶
In [3]:
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
In [4]:
mito_genes = np.array([name for name in adata.var_names if name.startswith('mt-')])
# for each cell compute fraction of counts in mito genes vs. all genes
adata.smp['percent_mito'] = np.sum(adata[:, mito_genes].X, axis=1).A1 / np.sum(adata.X, axis=1).A1
# add the total counts per cell as sample annotation to adata
adata.smp['n_counts'] = np.sum(adata.X, axis=1).A1
In [5]:
axs = sc.pl.scatter(adata, x='n_counts', y='percent_mito')
axs = sc.pl.scatter(adata, x='n_counts', y='n_genes')
In [6]:
print(adata.n_smps)
adata[adata.smp['n_genes'] < 2500, :]
adata = adata[adata.smp['percent_mito'] < 0.05, :]
print(adata.n_smps)
In [7]:
sc.write(data_folder+'pbmc3k_filtered_raw_log', sc.pp.log1p(adata, copy=True))
In [8]:
sc.pp.normalize_per_cell(adata, counts_per_cell_after=1e4)
filter_result = sc.pp.filter_genes_dispersion(
adata.X, min_mean=0.0125, max_mean=3, min_disp=0.5)
sc.pl.filter_genes_dispersion(filter_result)
In [9]:
adata = adata[:, filter_result.gene_subset]
sc.pp.log1p(adata)
In [10]:
adata_corrected = sc.pp.regress_out(adata, ['n_counts', 'percent_mito'], copy=True)
In [11]:
sc.pp.scale(adata_corrected, max_value=10)
In [12]:
sc.write(data_folder+'pbmc3k_corrected', adata_corrected)
Dimensionality Reduction¶
In [13]:
sc.tl.pca(adata_corrected)
ax = sc.pl.pca_scatter(adata_corrected, right_margin=0.2)
In [14]:
sc.pl.pca_variance_ratio(adata_corrected, log=0)
In [15]:
adata_corrected = sc.read(data_folder+'pbmc3k_corrected')
In [16]:
sc.tl.tsne(adata_corrected, n_pcs=15)
In [17]:
ax = sc.pl.tsne(adata_corrected)
Clustering¶
In [18]:
sc.tl.louvain(adata_corrected, flavor="igraph", recompute_graph=True, n_pcs=15)
In [19]:
plt.figure(figsize=[6,6])
ax = plt.subplot(111)
sc.pl.tsne(adata_corrected, color='louvain_groups', title="Louvain Clusters", ax=ax, show=0)
plt.show()
In [20]:
sc.write(data_folder+'pbmc3k_corrected', adata_corrected)
In [21]:
adata_corrected = sc.read(data_folder+'pbmc3k_corrected')
adata_raw = sc.read(data_folder+'pbmc3k_filtered_raw_log')
adata_raw.smp['louvain_groups'] = adata_corrected.smp['louvain_groups']
adata_raw.smp["X_tsne"] = adata_corrected.smp["X_tsne"]
sc.tl.rank_genes_groups(adata_raw, 'louvain_groups')
sc.pl.rank_genes_groups(adata_raw, n_genes=20)
In [22]:
diff_genes = sc.tl.rank_genes_groups(adata_raw, 'louvain_groups', copy=1, n_genes=400)
diff_genes_scores = pd.DataFrame(diff_genes.add["rank_genes_groups_gene_scores"])
diff_genes = pd.DataFrame(diff_genes.add["rank_genes_groups_gene_names"])
In [23]:
def top_clust_genes(clust_num, top_n=5):
print("Plotting Heatmaps for top "+str(top_n)+" genes in cluster #"+str(clust_num))
for gene_name in np.array(diff_genes.iloc[:top_n, clust_num]):
sc.pl.tsne(adata_raw, color=gene_name, title=gene_name+" Expression", color_map=cm1)
In [24]:
top_clust_genes(0, 10)
top_clust_genes(1, 10)
top_clust_genes(9, 10)
top_clust_genes(2)
top_clust_genes(7)
In [25]:
ax = sc.pl.violin(adata_raw, 'Ly6a', group_by='louvain_groups')
Gene Ontology Enrichment¶
In [26]:
import gzip
import Bio.UniProt.GOA as GOA
from goatools.base import download_go_basic_obo
from goatools import obo_parser, GOEnrichmentStudy
from goatools.godag_plot import plot_gos, plot_results, plot_goid2goobj
obo_fname = download_go_basic_obo()
go_db = obo_parser.GODag(obo_fname)
# File is a gunzip file, so we need to open it in this way
with gzip.open("goa_mouse.gaf.gz", 'rt') as arab_gaf_fp:
arab_funcs = {} # Initialise the dictionary of functions
# Iterate on each function using Bio.UniProt.GOA library.
for entry in GOA.gafiterator(arab_gaf_fp):
uniprot_id = entry.pop('DB_Object_ID')
arab_funcs[uniprot_id] = entry
gene_frame = pd.DataFrame(arab_funcs).T
In [27]:
def top_clust_ontologies(clust_num, top_n_ontologies=10, top_n_genes=200):
print("Plotting top "+str(top_n_ontologies)+" gene ontologies in cluster #"+str(clust_num)+" using the top "+str(top_n_genes)+" genes")
go_frame = pd.DataFrame([diff_genes.iloc[:top_n_genes, clust_num], diff_genes_scores.iloc[:top_n_genes, clust_num]]).T
go_frame.columns = ["GeneName", "GeneScore"]
go_frame = go_frame.merge(gene_frame[["DB_Object_Symbol", "GO_ID"]], how="left", left_on="GeneName", right_on="DB_Object_Symbol")
go_frame_sparse = go_frame.dropna()
missing_annot = len(go_frame.index) - len(go_frame_sparse.index)
if missing_annot > 0:
print("Missing annotations for "+str(missing_annot)+" of "+str(len(go_frame.index))+" genes")
go_frame_sparse = go_frame_sparse.assign(GO_Name= [go_db[x].name for x in go_frame_sparse["GO_ID"]])
go_scores = go_frame_sparse.groupby("GO_Name")["GeneScore"].sum().sort_values()[-top_n_ontologies:]
go_scores.plot.barh()
plt.title("Enriched Ontologies in Cluster #"+str(clust_num))
plt.show()
print("\n"+"-"*100+"\n")
return go_frame
In [28]:
_ = [top_clust_ontologies(x, 10, 100) for x in [0, 1, 9, 2, 7]]
In [29]:
pop = arab_funcs.keys()
assoc = {}
for x in arab_funcs:
name = arab_funcs[x]["DB_Object_Symbol"]
if name not in assoc:
assoc[name] = set()
assoc[name].add(str(arab_funcs[x]['GO_ID']))
pop = assoc.keys()
In [30]:
methods = ['fdr_bh']
g = GOEnrichmentStudy(pop, assoc, go_db,
propagate_counts=False,
alpha=0.05,
methods=methods)
In [31]:
def goatools_clust_ontologies(clust_num, top_n=200, pval=.01):
print("Plotting top "+str(top_n)+" gene ontologies in cluster #"+str(clust_num))
go_frame = pd.DataFrame([diff_genes.iloc[:, clust_num], diff_genes_scores.iloc[:, clust_num]]).T
go_frame.columns = ["GeneName", "GeneScore"]
go_frame = go_frame.merge(gene_frame[["DB_Object_Symbol", "GO_ID"]], how="left", left_on="GeneName", left_index=True, right_on="DB_Object_Symbol")
study_out = g.run_study(np.array(diff_genes.iloc[:top_n, clust_num]))
print([(r.p_fdr_bh, r.name, r.GO) for r in study_out if r.get_pvalue() < pval])
#g.print_summary(study_out, min_ratio=None, pval=pval)
sig_gene = [r for r in study_out if r.p_fdr_bh < pval]
sig_gene_GO = [r.GO for r in study_out if r.p_fdr_bh < pval]
#plot_results("write/"+data_folder+"clust_"+str(clust_num)+"__{NS}.png", sig_gene)
plot_gos("write/"+data_folder+"clust_"+str(clust_num)+".png",
sig_gene_GO, # Source GO ids
go_db,
goea_results=study_out, # use pvals for coloring
# We can further configure the plot...
study_items=6, # Only only 6 gene Symbols max on GO terms
items_p_line=3, # Print 3 genes per line
)
print("\n"+"-"*100+"\n")
return study_out
In [32]:
_ = [goatools_clust_ontologies(x, 400, .01) for x in [0, 1, 9, 2, 7]]